!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine hllc_x (fro, fmx, fmy, fmz, fbx, fby, fbz, fen,                                                                         &
   ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, po_l, ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, po_r, nx, ny, nz)
!==================================================================================================================================|

   use parameters
   implicit none

   integer(4), intent(in) :: nx, ny, nz
   real(8), dimension(nx, ny, nz), intent(in) :: ro_l, mx_l, my_l, mz_l, bx_l, by_l, bz_l, en_l, po_l
   real(8), dimension(nx, ny, nz), intent(in) :: ro_r, mx_r, my_r, mz_r, bx_r, by_r, bz_r, en_r, po_r
   real(8), dimension(nx, ny, nz), intent(inout) :: fro, fmx, fmy, fmz, fbx, fby, fbz, fen

   integer(4) :: ip1
   real(8) :: s_m, s_l, s_r
   real(8) :: ro_hll
   real(8) :: ro_1l, mx_1l, my_1l, mz_1l, bx_1l, by_1l, bz_1l, en_1l
   real(8) :: ro_1r, mx_1r, my_1r, mz_1r, bx_1r, by_1r, bz_1r, en_1r
   real(8) :: vx_l, vy_l, vz_l, pt_l, pr_l, v2_l, b2_l, s2_l, c2_l, cax2_l, cfx_l
   real(8) :: vx_r, vy_r, vz_r, pt_r, pr_r, v2_r, b2_r, s2_r, c2_r, cax2_r, cfx_r
   real(8) :: vx_1l, vy_1l, vz_1l, pt_1l
   real(8) :: vx_1r, vy_1r, vz_1r, pt_1r
   real(8) :: fro_l, fmx_l, fmy_l, fmz_l, fbx_l, fby_l, fbz_l, fen_l
   real(8) :: fro_r, fmx_r, fmy_r, fmz_r, fbx_r, fby_r, fbz_r, fen_r
   real(8) :: ex1, ex2, ex3

!----------------------------------------------------------------------------------------------------------------------------------|
   do k = 1, nz
   do j = 1, ny
   do i = 1, nx - 1
      ip1 = i + 1
!----------------------------------------------------------------------------------------------------------------------------------|
!     left vx, vy, vz, s, pt
!----------------------------------------------------------------------------------------------------------------------------------|
      vx_l = mx_l(i, j, k) / ro_l(i, j, k)
      vy_l = my_l(i, j, k) / ro_l(i, j, k)
      vz_l = mz_l(i, j, k) / ro_l(i, j, k)
      b2_l = bx_l(i, j, k) * bx_l(i, j, k) + by_l(i, j, k) * by_l(i, j, k) + bz_l(i, j, k) * bz_l(i, j, k)
      v2_l = vx_l * vx_l + vy_l * vy_l + vz_l * vz_l
      pr_l = (en_l(i, j, k) - v2_l * ro_l(i, j, k) / 2.0d0 - b2_l / 2.0d0) * gmm1
      c2_l = gm * pr_l / ro_l(i, j, k)
      s2_l = c2_l + b2_l / ro_l(i, j, k)
      cax2_l = bx_l(i, j, k) * bx_l(i, j, k) / ro_l(i, j, k)
      cfx_l = sqrt((s2_l + sqrt(s2_l * s2_l - 4.0d0 * c2_l * cax2_l)) / 2.0d0)

!----------------------------------------------------------------------------------------------------------------------------------|
!     right vx, vy, vz, s, pt
!----------------------------------------------------------------------------------------------------------------------------------|
      vx_r = mx_r(ip1, j, k) / ro_r(ip1, j, k)
      vy_r = my_r(ip1, j, k) / ro_r(ip1, j, k)
      vz_r = mz_r(ip1, j, k) / ro_r(ip1, j, k)
      b2_r = bx_r(ip1, j, k) * bx_r(ip1, j, k) + by_r(ip1, j, k) * by_r(ip1, j, k) + bz_r(ip1, j, k) * bz_r(ip1, j, k)
      v2_r = vx_r * vx_r + vy_r * vy_r + vz_r * vz_r
      pr_r = (en_r(ip1, j, k) - v2_r * ro_r(ip1, j, k) / 2.0d0 - b2_r / 2.0d0) * gmm1
      c2_r = gm * pr_r / ro_r(ip1, j, k)
      s2_r = c2_r + b2_r / ro_r(ip1, j, k)
      cax2_r = bx_r(ip1, j, k) * bx_r(ip1, j, k) / ro_r(ip1, j, k)
      cfx_r = sqrt((s2_r + sqrt(s2_r * s2_r - 4.0d0 * c2_r * cax2_r)) / 2.0d0)

!----------------------------------------------------------------------------------------------------------------------------------|
!     left and right function
!----------------------------------------------------------------------------------------------------------------------------------|
      fro_l = mx_l(i, j, k)
      fmx_l = vx_l * mx_l(i, j, k) + pr_l + b2_l / 2.0d0 - bx_l(i, j, k) * bx_l(i, j, k)
      fmy_l = vx_l * my_l(i, j, k) - bx_l(i, j, k) * by_l(i, j, k)
      fmz_l = vx_l * mz_l(i, j, k) - bx_l(i, j, k) * bz_l(i, j, k)
      fbx_l = po_l(i, j, k)
      fby_l = vx_l * by_l(i, j, k) - bx_l(i, j, k) * vy_l
      fbz_l = vx_l * bz_l(i, j, k) - bx_l(i, j, k) * vz_l
      fen_l = (en_l(i, j, k) + pr_l + b2_l / 2.0d0) * vx_l -                                                                       &
         (bx_l(i, j, k) * vx_l + by_l(i, j, k) * vy_l + bz_l(i, j, k) * vz_l) * bx_l(i, j, k)

      fro_r = mx_r(ip1, j, k)
      fmx_r = vx_r * mx_r(ip1, j, k) + pr_r + b2_r / 2.0d0 - bx_r(ip1, j, k) * bx_r(ip1, j, k)
      fmy_r = vx_r * my_r(ip1, j, k) - bx_r(ip1, j, k) * by_r(ip1, j, k)
      fmz_r = vx_r * mz_r(ip1, j, k) - bx_r(ip1, j, k) * bz_r(ip1, j, k)
      fbx_r = po_r(ip1, j, k)
      fby_r = vx_r * by_r(ip1, j, k) - bx_r(ip1, j, k) * vy_r
      fbz_r = vx_r * bz_r(ip1, j, k) - bx_r(ip1, j, k) * vz_r
      fen_r = (en_r(ip1, j, k) + pr_r + b2_r / 2.0d0) * vx_r -                                                                     &
         (bx_r(ip1, j, k) * vx_r + by_r(ip1, j, k) * vy_r + bz_r(ip1, j, k) * vz_r) * bx_r(ip1, j, k)

!----------------------------------------------------------------------------------------------------------------------------------|
!     left and right var*
!----------------------------------------------------------------------------------------------------------------------------------|
      s_l = min (vx_l - cfx_l, vx_r - cfx_r)
      s_r = max (vx_l + cfx_l, vx_r + cfx_r)
      pt_l = pr_l + b2_l / 2.0d0
      pt_r = pr_r + b2_r / 2.0d0
      ex1 = (s_r - vx_r) * ro_r(ip1, j, k) - (s_l - vx_l) * ro_l(i, j, k)
      if (ex1 .eq. 0.0d0) ex1 = eps
      s_m = ((s_r - vx_r) * ro_r(ip1, j, k) * vx_r - (s_l - vx_l) * ro_l(i, j, k) * vx_l -                                         &
         pt_r + pt_l - bx_l(i, j, k) * bx_l(i, j, k) + bx_r(ip1, j, k) * bx_r(ip1, j, k)) / ex1

      ex3 = s_l - s_m
      if (ex3 .eq. 0.0d0) ex3 = eps
      ex1 = (s_l - vx_l) / ex3
      ex2 = s_r - s_l
      if (ex2 .eq. 0.0d0) ex2 = eps
      ro_1l = ro_l(i, j, k) * ex1
      ro_hll = (s_r * ro_r(ip1, j, k) - s_l * ro_l(i, j, k) - (fro_r - fro_l)) / ex2
      if (ro_hll .eq. 0.0d0) ro_hll = eps
      bx_1l = (s_r * bx_r(ip1, j, k) - s_l * bx_l(i, j, k) - (fbx_r - fbx_l)) / ex2
      by_1l = (s_r * by_r(ip1, j, k) - s_l * by_l(i, j, k) - (fby_r - fby_l)) / ex2
      bz_1l = (s_r * bz_r(ip1, j, k) - s_l * bz_l(i, j, k) - (fbz_r - fbz_l)) / ex2
      vx_1l = (s_r * mx_r(ip1, j, k) - s_l * mx_l(i, j, k) - (fmx_r - fmx_l)) / ex2 / ro_hll
      vy_1l = (s_r * my_r(ip1, j, k) - s_l * my_l(i, j, k) - (fmy_r - fmy_l)) / ex2 / ro_hll
      vz_1l = (s_r * mz_r(ip1, j, k) - s_l * mz_l(i, j, k) - (fmz_r - fmz_l)) / ex2 / ro_hll
      mx_1l = s_m * ro_1l
      my_1l = my_l(i, j, k) * ex1 - (bx_1l * by_1l - bx_l(i, j, k) * by_l(i, j, k)) / ex3
      mz_1l = mz_l(i, j, k) * ex1 - (bx_1l * bz_1l - bx_l(i, j, k) * bz_l(i, j, k)) / ex3
      pt_1l = ro_l(i, j, k) * (s_l - vx_l) * (s_m - vx_l) + pt_l - bx_l(i, j, k) * bx_l(i, j, k) + bx_1l * bx_1l
      en_1l = en_l(i, j, k) * ex1 + (pt_1l * s_m - pt_l * vx_l -                                                                   &
         (bx_1l * (vx_1l * bx_1l + vy_1l * by_1l + vz_1l * bz_1l) -                                                                &
          bx_l(i, j, k) * (vx_l * bx_l(i, j, k) + vy_l * by_l(i, j, k) + vz_l * bz_l(i, j, k)))) / ex3

      ex3 = s_r - s_m
      if (ex3 .eq. 0.0d0) ex3 = eps
      ex1 = (s_r - vx_r) / ex3
      ro_1r = ro_r(ip1, j, k) * ex1
      bx_1r = bx_1l
      by_1r = by_1l
      bz_1r = bz_1l
      vx_1r = vx_1l
      vy_1r = vy_1l
      vz_1r = vz_1l
      mx_1r = s_m * ro_1r
      my_1r = my_r(ip1, j, k) * ex1 - (bx_1r * by_1r - bx_r(ip1, j, k) * by_r(ip1, j, k)) / ex3
      mz_1r = mz_r(ip1, j, k) * ex1 - (bx_1r * bz_1r - bx_r(ip1, j, k) * bz_r(ip1, j, k)) / ex3
      pt_1r = ro_r(ip1, j, k) * (s_r - vx_r) * (s_m - vx_r) + pt_r - bx_r(ip1, j, k) * bx_r(ip1, j, k) + bx_1r * bx_1r
      en_1r = en_r(ip1, j, k) * ex1 + (pt_1r * s_m - pt_r * vx_r -                                                                 &
         (bx_1r * (vx_1r * bx_1r + vy_1r * by_1r + vz_1r * bz_1r) -                                                                &
          bx_r(ip1, j, k) * (vx_r * bx_r(ip1, j, k) + vy_r * by_r(ip1, j, k) + vz_r * bz_r(ip1, j, k)))) / ex3

!----------------------------------------------------------------------------------------------------------------------------------|
!     flux
!----------------------------------------------------------------------------------------------------------------------------------|
      if (s_l .gt. 0.0d0) then
         fro(i, j, k) = fro_l
         fmx(i, j, k) = fmx_l
         fmy(i, j, k) = fmy_l
         fmz(i, j, k) = fmz_l
         fbx(i, j, k) = fbx_l
         fby(i, j, k) = fby_l
         fbz(i, j, k) = fbz_l
         fen(i, j, k) = fen_l
      elseif (s_l .le. 0.0d0 .and. s_m .gt. 0.0d0) then
         fro(i, j, k) = fro_l + s_l * (ro_1l - ro_l(i, j, k))
         fmx(i, j, k) = fmx_l + s_l * (mx_1l - mx_l(i, j, k))
         fmy(i, j, k) = fmy_l + s_l * (my_1l - my_l(i, j, k))
         fmz(i, j, k) = fmz_l + s_l * (mz_1l - mz_l(i, j, k))
         fbx(i, j, k) = fbx_l + s_l * (bx_1l - bx_l(i, j, k))
         fby(i, j, k) = fby_l + s_l * (by_1l - by_l(i, j, k))
         fbz(i, j, k) = fbz_l + s_l * (bz_1l - bz_l(i, j, k))
         fen(i, j, k) = fen_l + s_l * (en_1l - en_l(i, j, k))
      elseif (s_m .le. 0.0d0 .and. s_r .gt. 0.0d0) then
         fro(i, j, k) = fro_r + s_r * (ro_1r - ro_r(ip1, j, k))
         fmx(i, j, k) = fmx_r + s_r * (mx_1r - mx_r(ip1, j, k))
         fmy(i, j, k) = fmy_r + s_r * (my_1r - my_r(ip1, j, k))
         fmz(i, j, k) = fmz_r + s_r * (mz_1r - mz_r(ip1, j, k))
         fbx(i, j, k) = fbx_r + s_r * (bx_1r - bx_r(ip1, j, k))
         fby(i, j, k) = fby_r + s_r * (by_1r - by_r(ip1, j, k))
         fbz(i, j, k) = fbz_r + s_r * (bz_1r - bz_r(ip1, j, k))
         fen(i, j, k) = fen_r + s_r * (en_1r - en_r(ip1, j, k))
      elseif (s_r .le. 0.0d0) then
         fro(i, j, k) = fro_r
         fmx(i, j, k) = fmx_r
         fmy(i, j, k) = fmy_r
         fmz(i, j, k) = fmz_r
         fbx(i, j, k) = fbx_r
         fby(i, j, k) = fby_r
         fbz(i, j, k) = fbz_r
         fen(i, j, k) = fen_r
      endif
   enddo
   enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine hllc_x
